% num = match3(image1, image2,image3)
%
% This function takes two images, finds their SIFT features, and
%   displays lines connecting the matched keypoints.  A match is accepted
%   only if its distance is less than distRatio times the distance to the
%   second closest match.
% It returns the number of matches displayed.

function [num loc1 loc2 loc3 des1 des2 des3 matched1 matched2] = match3(image1, image2, image3)

% Find SIFT keypoints for each image
[im1, des1, loc1] = sift(image1);
[im2, des2, loc2] = sift(image2);
[im3, des3, loc3] = sift(image3);

% For efficiency in Matlab, it is cheaper to compute dot products between
%  unit vectors rather than Euclidean distances.  Note that the ratio of 
%  angles (acos of dot products of unit vectors) is a close approximation
%  to the ratio of Euclidean distances for small angles.
%
% distRatio: Only keep matches in which the ratio of vector angles from the
%   nearest to second nearest neighbor is less than distRatio.
distRatio = 0.6;   

% For each descriptor in the first image, select its match to second image.
des2t = des2';                          % Precompute matrix transpose
des3t = des3';                          % Precompute matrix transpose
matched1 = zeros(size(des1,1),1);
for i = 1 : size(des1,1)
   dotprods1 = des1(i,:) * des2t;        % Computes vector of dot products
   [vals1,indx1] = sort(acos(dotprods1));  % Take inverse cosine and sort results
   % Check if nearest neighbor has angle less than distRatio times 2nd.
	if (vals1(1) < distRatio * vals1(2))
        matched1(i) = indx1(1);
	end
end

matched2 = zeros(size(des2,1),1);
for i= 1:size(des2,1)
    dotprods2 = des2(i,:) * des3t;        % Computes vector of dot products
    [vals2,indx2] = sort(acos(dotprods2));  % Take inverse cosine and sort results
    if (vals2(1) < distRatio * vals2(2))
        matched2(i) = indx2(1);       
    end
end

% matched3 = zeros(size(des1,1),1);
% for i = 1:size(des1,1)
%     for j=1:size(des2,1)
%         %matched
%         if matched1(i) == matched2(j) && matched1(i)~=0
%             matched3(i)=j;
%         end
%     end
% end

% Create a new image showing the two images side by side.
im4 = appendimages(appendimages(im1,im2),im3);

% Show a figure with lines joining the accepted matches.
figure;
imshow(im4);
hold on;
cols1 = size(im1,2);
num = 0;
for i = 1: size(des1,1)
    if (matched1(i) > 0)
        if (matched2(matched1(i)) > 0)
            %line([loc1(i,2) loc2(match(i),2)+cols1], ...
            %     [loc1(i,1) loc2(match(i),1)], 'Color', 'c');
            plot([loc1(i,2) loc2(matched1(i),2)+cols1 loc3(matched2(matched1(i)),2)+cols1*2], ...
                [loc1(i,1) loc2(matched1(i),1) loc3(matched2(matched1(i)),1)], '.b');
            num = num + 1;
        end
    end
end
hold off;
fprintf('Found %d matches.\n', num);




